Lattice dynamics of anharmonic solids from first principles 



O 
(N 

o 
O 



o 



> 
o 
o\ 
in 

rn 
O 



X 



O. Hellman,^ I. A. Abrikosov/ and S. I. Siniak^ 

^Department of Physics, Chemistry and Biology (IFM), 
Linkoping University, SE-581 83, Linkoping, Sweden. 

An accurate and easily extendable method to deal with lattice dynamics of solids is offered. It is 
based on first-principles molecular dynamics simulations and provides a consistent way to extract 
the best possible harmonic - or higher order - potential energy surface at finite temperatures. 
It is designed to work even for strongly anharmonic systems where the traditional quasiharmonic 
approximation fails. The accuracy and convergence of the method are controlled in a straightforward 
way. Excellent agreement of the calculated phonon dispersion relations at finite temperature with 
experimental results for bcc Li and bcc Zr is demonstrated. 



Ability to predict phase equilibria and structural trans- 
formations in solids under pressure and temperature is of 
vital importance to both science and industry. The den- 
sity functional theory (DFT) 1] allows one to calculate 
thermodynamic properties of solids from first principles, 
i.e. without any adjustable parameters or information 
from experiment. Though it is a ground-state theory the 
crucial effect of thermal vibrations of atoms can inprinci- 
ple be included via the harmonic approximation[2|. The 
DFT-based calculation of corresponding phonon disper- 
sion relations is nowadays routine^] . The widely used so- 
called quasiharmonic approximation Q for the free energy 
employing zero temperature DFT calculations of volume- 
dependent phonon frequencies is well developed ,4]. How- 
ever, it assumes that the potential energy surface does 
not depend on temperature and is limited to systems 
where anharmonic effects can be neglected. In particular 
it is known to dramatically fail when applied to crystal 
structures which are dynamically unstable at K, i.e. 
their phonon spectra contain imaginary frequencies, but 
which stabilized due to the anharmonic contributions to 
the free energy at elevated temperatures [5|. 

The failure of the harmonic theory mostly arises from 
the K Taylor expansion of the potential energy sur- 
face up to just the quadratic term, which at temper- 
atures close to the melting point becomes increasingly 
inaccurate. Consequently, theoretical treatment of an- 
harmonic effects has attracted great interest [6|. The self- 
consistent phonon approach was formulated by Born|7|l, 
and developed by Born and Hooton[8, 9]. Choquard lOj 
and Plakida[llj addressed the problem employing the 
perturbation theory based on the so-called double-time 
Green functions. However, such methods are difficult 
to combine with DFT-based first-principles calculations. 
The expansion to higher-order terms within the standard 
techniques is cumbersome and currently unfeasible, in 
particular due to issues of numerical precision. 



On the other hand, recently Souvatzis et al.[l2| have 
introduced the so-called self-consistent ab initio lattice 
dynamics approach, which is conceptually similar to the 
so-called renormalized harmonic approximation The 
theory introduces anharmonic effects through the tem- 



perature dependence of phonon frequencies mimicking in- 
teraction of phonons. The price is the additional approx- 
imations, such as fixed amplitudes of the vibrations and 



a limited phase space of allowed excitations. In Ref. [12 1 
experimental dispersion relations for the group IV met- 
als have in general been reproduced, though substantial 
error has been reported for a determination of transition 
temperatures between different phases [l3| . 

Fourier analyzing the velocity autocorrelation function 
from molecular dynamics (MD) represents an alterna- 
tive approach to determine the phonon dispersion rela- 
tions at finite temperature. It works for classical MD 
simulations (T3|. but the system sizes currently accessible 
to ah initio MD (AIMD) makes this approach prone to 
finite size effects and can thus not be used with meaning- 
ful accuracy. The common practical way of dealing with 
anharmonic effects in DFT based simulations of phase 
transformations is combine AIMD with a thermodynamic 
integration using the harmonic energy as the starting 
point [Til. [l6|. In the standard harmonic approximation, 
however, the imaginary phonon frequencies mean that 
the free energy is not defined and therefore the starting 
point is missing. Therefore often an artificial model with 
a fitted classical potential is used instead [l5|. 

We propose a method that intrinsically solves all the 
issues mentioned above. It uses ah initio molecular dy- 
namics simulations and provides a consistent and compu- 
tationally easy way to extract the best possible harmonic 
- or higher order - potential energy surface at finite tem- 
peratures to study the lattice dynamics and thermody- 
namic properties of the system under consideration. At 
present we show that already the second order terms are 
sufficient to accurately describe systems where the anhar- 
monic effects are known to be prevalent and as examples 
we present calculations for the strongly anharmonic bcc 
Li and Zr. We stress, however, that if the second-order 
approach proves insufficient for some complex system it 
is easily extended to any order. Further, our method al- 
lows one to control the accuracy and convergence in a 
straightforward way, in contrast to most of the previous 
methods. 

For clarity the equations are derived for a monoatomic 
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system but it is easy to extend them to the general case. 
Usually the harmonic approximation is given as a Tay- 
lor expansion of the potential energy surface (at K) 
in atomic displacements u, truncating after the second 
order terms: 



(1) 

where the second derivatives of the energy surface U with 
respect to displacements u are called the force constant 
matrices, D'^^j^. Here Hi and Rj are ideal atomic posi- 
tions and n and z/ correspond to the Cartesian coordi- 
nates. The established way to find the force constant 
matrices D is through the linear response 3j or the small 
displacement mcthodT7|. In these methods the crystal is 
considered in its ground state and displacements of indi- 
vidual atoms, either via perturbation theory or a super- 
cell approach, are introduced. In this work however, we 
want to obtain the best possible harmonic potential for a 
given temperature and volume, derived from a thermally 
excited crystal. 

In an expansion in the second order the forces are given 
by 



(2) 



We want to find the effective force constant matrix D 
that best represent the real forces. To do this we run Nt 
time steps of ab initio molecular dynamics simulations 
for a supercell containing Na atoms, and at each time 
step t we store forces and displacements. We then seek 
to minimise 
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(3) 



where F* is the force acting on atom i at time step t 
from the simulation, and F ■ is the force given by Eq. [2] 
with displacements taken from time step t. If the number 
of time steps are larger than 3Na we have an overdeter- 
mined system of equations for D. We can then obtain 
the linear least squares solution of D. Having converged 
with respect to the number of time steps we have defined 
a force constant matrix at temperature T we can obtain 
the phonon dispersion relations and free energy. 

To monitor the quality of the force constant matrix for 
bcc Li we determine it from successively longer simula- 
tions, and show the result in Fig. [1] It clearly demon- 
strates that the suggested scheme allows us to derive a 
unique and well converged effective force constant matrix 
D. 

We explicitly note that the forces, displacements and 
force constant matrix in this example are obtained at 
temperature T = 300K, whereas the quasiharmonic ap- 
proximation would derive Z? at K and from that obtain 
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FIG. 1. (color online) Difference between the free energy 
(solid line, left axis) and ||D^^(Rj)|| (dashed line, right 
axis), the Frobenius norm of the nearest neighbour interac- 
tion derived from increasingly longer simulations for bcc Li 
at r = 300K. As a reference the results obtained from simu- 
lations with 25000 time steps are used, which are considered 
to give the converged results (0.43 eV/A^ for ||_D"(Ri)||). 



the free energy as a function of temperature. Note that 
in our method the force constant matrix D^^ is defined 
per pair ij in the supercell. However, a careful remap- 
ping allows one to define it for each lattice vector, that 
is 



(4) 



where a, (5 are indices to atoms in a suitable unit cell of 
choice (smaller or equal to the supercell), R^ is vector i in 
coordination shell k that contains only pairs of type a, /3. 
To increase numerical stability we use the point group of 
the coordination shells as well as symmetry properties 
of the force constant matrices. All the symmetrization 
procedures are similar to those detailed in Alfe et al. , 
although used for a different purpose. The small dis- 
placement method uses as few displacements as possible 
and then constructs all force constant matrices through 
symmetry relations, putting high accuracy demands on 
the calculated forces. Our method yields all force con- 
stant matrices that are then forced to obey the symmetry 
relations of the lattice, reducing the need for numerical 
accuracy in the calculated forces. 

16[ we analyzed the forces to 



In the spirit of Ref. 



asses the errors in our method. We analysed forces from 
a subset of configurations, calculated with different level 
of accuracy. F^, forces with an accuracy typical for first 
principles MD simulations were compared to fully 
converged results[3| ■ We found that the norm ratio, 
oscillated around unity (« 1.0015), with 
standard deviation of 0.007. This, as well as the results 
from several test calculations and data presented in Fig. 
[U allows us to draw the conclusion that since the error 
in forces is normally distributed around the "true" value, 
our method of least square fitting of many configurations 
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appearing during the MD simulations can reproduce force 
constant matrices with high accuracy. 

The fact that we do not do a Taylor expansion around 
equilibrium at K but produce the potential energy sur- 
face according to the least square fit of a quadratic form 
at a particular finite temperature allows one to apply 
our method to systems where the quasiharmonic approx- 
imation traditionally fails. We notice that the standard 
quasiharmonic approximation has a temperature range 
at low temperatures where it works well for dynamically 
stable systems. Our method moves this window to any 
temperature of interest. Although not explicitly anhar- 
monic, it gives us truly the best harmonic fit to the fully 
anharmonic energy landscape and therefore contains the 
information about anharmonism implicitly through the 
temperature dependent vibrational frequencies. 

At this point it should be stressed once again that 
this method allows for an easy, though a more time- 
consuming, extension to explicitly handle anharmonism 
by including more terms in the expansion in Eq. ([1}. 

We chose bcc Li and Zr to demonstrate the advantages 
of the proposed technique. Li, long thought of as a sim- 
ple metal, has a complex phase diagram. It undergoes a 
phase transition from the bcc structure at room temper- 
ature to the close-packed R9 structure below 70 K[19|. 
With pressure it behaves anomalo usly transforming to 
the low symmetry Cmca structure 20| . Quasiharmonic 
phonon dispersion relations for the bcc structure show 
that it is dynamically unstable at K, and the free en- 
ergy is not defined. Zr has long been used as a model 
system for martensitic phase transitions and is a well 
known example of a strongly anharmonic solid. The na- 
ture and origin of the stabilization of the bcc phase has 
been discussed in numerous previous studies [2 

All electronic structure calculations were carried out 
with the projector-augmented wave (PAW) method as 
implemented in the code VASP|22h25|. We used a 128 
atom bcc supercell (4x4x4) for the MD simulations. 
For the BZ integration we used a 2 x 2 x 2 k-point mesh 
and Fermi smearing corresponding to the simulation tem- 
perature. Exchange and correlation effects were treated 
using the Perdew-Burke-Ernzerhof|26 functional form. 
A plane wave cutoff of 140 eV was used for Li and 154 eV 
for Zr. Both systems were considered at their theoretical 
equilibrium lattice parameter. 

We used a 2 fs time step, which is suitable for both 
systems, and set the temperature with the Nose-Hoover 
thermostat 271 128| . To fully check the convergence of our 
method we ran the calculations for a total of 44000 time 
steps (88 ps) after the initial equilibration and extracted 
the force constant matrix at fixed time intervals. After 
about 40 ps the free energy was converged to below 0.5 
me V/ atom (see Fig. [T]), which is an accuracy exceeding 
that of the underlying DFT approximations. 

The success of the suggested method on the quantita- 
tive level can be judged from a direct comparison of the 



calculated phonon dispersion relations (easily extracted 
from the calculated force constant matrix) to experimen- 
tal. In Fig [2] calculated dispersion relations at OK, room 
temperature and experimental results are shown for bcc 
Li. The K dispersion relations, obtained from the 
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FIG. 2. (color online) Phonon dispersion relations in bcc Li 
along high-symmetry directions. The symbols are experimen- 
tal values (293 K)[23|, the solid lines correspond to calcula- 
tions carried out at 300 K with the method proposed in this 
work and the dashed black hues to the K harmonic calcu- 
lations. 

quasiharmonic approximation reveal imaginary frequen- 
cies along the T — N direction. Using our method at finite 
temperature all imaginary frequencies disappear and we 
have an excellent agreement with experimental values. 

Looking at Zr in Fig. [3] we see once again excellent 
agreement of the results obtained by our method with 
experimental values. Worth noting is the near perfect 
location of the soft mode at the so-called cj-point in the 
H-P direction. It is indicated by the red vertical line. 
This mode is of crucial importance for the bcc to lo phase 
transition in Zr and it has been difficult to reproduce in 
previous calculations [l^. 

The phonon dispersion relations, while interesting, are 
not the main goal of this work. We want a solid method 
to deal with lattice dynamics for strongly anharmonic 
dynamically unstable systems. To test this we calculated 
the Gibbs free energy surface for bcc and hep Zr, and in 
Fig. |3]we present the calculated bcc- hep phase diagram. 
The dynamically unstable bcc free energies were calcu- 
lated from the phonon dispersion relations as detailed 
above, on a grid of six volumes and five temperatures 
(volume ±20%, temperature 100-1700K). 

For dynamically stable hep phase we used the quasi- 
harmonic approximation jT3| with phonon dispersion re- 
lations obtained on six volumes. We observe excellent 
agreement with experimentally determined line for the 
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FIG. 3. (color online) Phonon dispersion relations for bcc Zr. 
Solid lines correspond to calculations at 1300 K, dashed lines 
to the quasiharmonic results and symbols to experimental val- 
ues from Helming et al.[30l|fcircles) and Stassis et al.[3Hffilled 
circles) . The dotted vertical line is at q — (|,|,|) The ob- 
served softening at this point, experimental as well as theo- 
retical, is important for the bcc-cj phase transition. 



hep to bcc transition over the whole interval of pressures 
and temperatures. 

In summary we have developed a thorough, accurate 
and easily extendable formalism to deal with lattice dy- 
namics at high temperatures where the traditional quasi- 
harmonic approximation fails. Excellent agreement with 
experimental results for bcc Li and bcc Zr indicates the 
usefulness of the method. 
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